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We have investigated the "weak chaos" exponent to see if it can be considered as a classification 
parameter of different sandpile models. Simulation results show that "weak chaos" exponent may 
be one of the characteristic exponents of the attractor of deterministic models. We have shown 
that the (Abelian) BTW sandpile model and the (non Abelian) Zhang model posses different "weak 
chaos" exponents, so they may belong to different universality classes. We have also shown that 
stochasticity destroys "weak chaos" exponents' effectiveness so it slows down the divergence of 
nearby configurations. Finally we show that getting off the critical point destroys this behavior of 
deterministic models. 
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I. INTRODUCTION 

Bak, Tang, and Wissenfeld (BTW) introduced the con- 
cept of self-organized criticality (SOC) and the so-called 
BTW sandpile model as a description of power spacial 
and temporal correlations observed in a wide range of 
natural phenomenon^. During the past two decades, 
more sandpile models were introduced by different vari- 
ations of the main paradigm of SOC, the BTW model, 
in order to gain more realistic models. These models 
differ in some properties such as: discrete or continuous 
height variable, Abelian or non Abelian toppling rule, 
stochastic or deterministic toppling rule, directed or non 
directed (or even non directed on average) toppling cur- 
rent, sticky or non-sticky grains, and etc. We now have a 
large number of different sandpile- like models each model 
having it's own set of critical exponents. The number of 
these models are exceeding but after about two decades, 
a little is known about their universal classification. Al- 
though some numerical studies are done to address the 
universality of the critical behavior of different sandpile 
models but these studies are in contradiction with each 
other. For example Manna's classification^ is in contrast 
with Bcn-Hur and Biham's^ one, and the latter is in con- 
tradiction with Chessa et al£- classification. Also more 
recently fixed energy sandpiles (FES) are introduced in 
order to study the critical behavior of SOC and classi- 
fying sandpile models which has not gained any serious 
success yet and even it does not seem to be a successful 
career (for example see£). 

The first classification was done by Manna. He put his 
model and the BTW model in the same universality class 
ascribing the observed difference in critical exponents to 
finite size effects^. This result was verified by Grass- 
berger and Manna^. Another effort was done by Daz- 
Guilera and Corral based on renormalization group which 
resulted in classifying BTW and Zhang model in the 
same clas o 9 ' 10 (which confirmed Zhang's conjecture^). 
Bcn-Hur and Biham4 studied the most complete set of 
critical exponents {s,a,t,d,r,p}, based on the evolu- 
tion of conditional expectation values (See Christensen et 
al*£). These exponents are related to: (s) the size of an 
avalanche, (a) the area of an avalanche, (t) the time du- 
ration of an avalanche, (d) the maximal distance between 



the origin and the sites that an avalanche cluster touches, 
(r) the radius of gyration of an avalanche cluster, and (p) 
the perimeter of an avalanche cluster. They showed that 
sandpile models are classified in three groups of non di- 
rected, non directed on average, and directed models^. 
As a result of their study, BTW and Zhang models be- 
long to the same universality class (non directed) which 
Manna model (as a non directed on average model) does 
not and directed models belong to another class. Chessa 
et al. made some systematic corrections on Bcn-Hur and 
Biham's method and claimed that both stochastic and 
deterministic sandpile models belong to the same univer- 
sality class^. 



On the other hand, Bak and Chen had investigated 
the chaotic behavior of a block-spring model (which was 
introduced for simulating earthquake dynamics)^. They 
have shown that although the largest Lyapanuv exponent 
of this model is zero but nearby configurations separate 
in a power- law manner and they called it "weak chaos"—. 
Based on this study, Bak, Tang, and Chen (BTC) conjec- 
tured that SOC takes the system to the border of chaos. 
They also argued that this behavior is not because of ex- 
ponential sensitivity to initial conditions but the critical 
fluctuations of the system. They conjectured that in this 
manner, "weak chaos" is another aspect of the critical- 
ity of the attractor so the "weak chaos" exponent is a 
characteristic exponent of the syste m 13 ' . 



Although Vicira and Lichtcnberg found a counter 
example for BTC's conjecture^, we have numerically 
checked it for different sandpile models. We have shown 
that BTC's conjecture is truly verified in BTW (also in 
CBTW) and Zhang models but it's not true in Manna 
stochastic model. Dhar-Ramaswamy directed model 
behaves more complicated with at least two different 
regimes. 



In this paper we first define some sandpile models and 
discuss the time evolution of nearby points in their config- 
uration space. We will finally discuss the same behavior 
in the off critical regime. 
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FIG. 1: Evolution of H(i) for BTW model on three differ- 
ent lattices on Log scale for L — 256 and 10000 samples. 
The down-right, up-left, and the main graphs correspond to 
square, honeycomb, and triangular lattices correspondingly. 
Lattice effects which slows down H(t) in the first few time 
steps is mostly seen on square lattice. 



FIG. 2: Evolution of H(t) for Zhang model on three differ- 
ent lattices on Log scale for L = 256 and 10000 samples. 
The main, down-right, and the up-left graphs correspond to 
square, triangular, and honeycomb lattices correspondingly. 
Lattice effects which slows down H(t) in the first few time 
steps is mostly seen on honeycomb lattice. 



II. DEFINING THE MODELS 

We have studied sandpile models on two dimensional 
square, triangular, and honeycomb lattices. A height 
variable is assigned to each site (i,j) of the lattice which 
can be discrete or continuous depending on the model. 
This height variable could be interpreted as "energy". 
At each time step, sand is added to a randomly chosen 
site. Whenever the height of a site Zij, exceeds the crit- 
ical height (zij > z c ), the site would relax through the 
related toppling rule. Relaxation of a site would cause 
other sites to become unstable so they would topple and 
a chain reaction called an avalanche continues until all 
sites become stable. The rate of energy injection is so 
slow that an unstable configuration will relax before next 
grain is added. 

BTW model: the height variable of this model is dis- 
crete and its critical height is 4, 6, and 3 on square, tri- 
angular, and honeycomb lattices correspondingly. When 
Zi,j > z c the (i,j)th site relaxes: 
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which nn means the nearest neighbors. The boundary 
condition of this model is chosen to be open on all bound- 
aries. 

Zhang model: the height variable in this model is con- 
tinuous and its critical height can always be taken to be 1. 
At each time step a random amount of sand 5, which we 
take it to be in the set (0, 0.25), is added to a randomly 
chosen site. If z 2J > 1 this site topples by: 



nnn is the number of nearest neighbors 
ditions in this model is also open on all 
Continuous BTW model (CBTW) 

continuous version of BTW modeli&. Aj 
variable is continuous, the critical heig 
to be 1. At each time step an amount 
< 8 < 0.25, is added to a randomly 
the toppling rule is given via 
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nnn is again the number of nearest neighbors. 
Manna model: in this stochastic model the critical 
height is taken to be 2 on square and triangular lattices. 
On the square lattice when an unstable site is going to 
topple, it cither gives its left and right or its up and down 
neighbors one unit of sand each with equal probability. 
On the triangular lattice when a site topples, it gives one 
grain of sand to each facing neighbors with a probability 
of 1/3. The boundary conditions on this model is chosen 
to be open on all boundaries. 

Directed model: Dhar-Ramaswamy modelil is referred 
as the directed model. This model is defined on a square 
lattice in the (1,1) direction. The critical height is 2 and 
when a site topples it gives one sand to each two down 
neighbors. The vertical boundary condition is cylindri- 
cal, sands arc added to the system from the most top row 
and leave it from the lowest row. 



III. CHAOS IN SANDPILE MODELS 

To study chaos in different sandpile models we have 
monitored the time evolution of the distance of two 
nearby configurations. The Hamming distance of two 





100 120 140 




time 



time 



FIG. 3: Evolution of H(t) for Manna model on square and 
triangular lattices on Log scale for L — 256 and 5000 samples. 
The main and up-left graphs correspond to square and trian- 
gular lattices correspondingly. The down-right graph shows a 
logarithmic fit for this model on square lattice. 



FIG. 4: Evolution of Hamming distance for Dhar- 
Ramaswamy directed model for L — 128 and 2000 samples. 
This model shows a more complicated behavior. There are 
two power-law regimes with 0.64 ± 0.04 and 0.15 ± 0.03 cor- 
respondingly 



configurations at time t which is defined by 



H(t) = {J2hAt)-4At)?Y 



(4) 



is used as the distance of two sandpiles^. We have first 
provided a configuration Z which its height variables 
are denoted by {zij}, then we manipulate the heights 
of about 0.0005 of the whole sites chosen randomly, to 
prepare a nearby configuration Z' with heights {z[ ■}. In 
each time step an amount of sand is added to a similar 
randomly chosen site of Z and Z' and H(t) is calculated 
after relaxation of both sandpiles. 

BTW model: figure [T] shows the evolution of H(t) for 
this model on square, triangular, and honeycomb lattices. 
Comparing H(t) on different lattices shows that (ignor- 
ing a few first time steps which is attributed to lattice 
effects) 



#(i)~^, 
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With (3 sgU are = 0.74 ±0.03, ^triangular = 0.76 ±0.03, and 
Phoneycomb = 0.76 ± 0.03. 

Simulations of this model is done for L = 
64, 100, 128, 200, 256, 300, 400, 512, 1024 and for L > 200 
the exponents are independent of L (for smaller sandpiles 
the exponents increase slightly by system size). As it's 
seen in fig [I] and table HI the "weak chaos" exponent does 
not depend on lattice geometry. This is a good evidence 
for for considering it as universal property of this model. 

The behavior of CBTW model is exactly the same as 
this model as expected. 

Zhang model: figure [2] shows the evolution of H(t) for 
this model again on square, triangular, and honeycomb 
lattices. Comparing these three shows that (ignoring the 
first few points) Zhang model obeys BTC's conjecture 

With exponents (3 sqU are = 0.96±0.05, ^triangular = 0.93± 



Model Lattice type 


size 


P 


BTW 


square 


256 


0.74 ±0.03 


BTW 


triangular 


256 


0.76 ±0.03 


BTW 


honeycomb 


256 


0.76 ± 0.03 


Zhang 


square 


256 


0.96 ±0.05 


Zhang 


triangular 


256 


0.93 ±0.05 


Zhang 


honeycomb 


256 


0.94 ± 0.05 



TABLE I: Weak chaos exponent for BTW and Zhang mod- 
els on square, triangular, and honeycomb lattices. Since lattice 
effects can not be separated from power-law regime in a stan- 
dard way, the errors are reported on the base of simulating 
each model five times each time containing 5000 samples. 



0.05, and Phoneycomb = 0.94 ±0.05. Again the exponents 
are not lattice dependent. 

Manna model: as it's seen in figure [3] the evolution of 
H{t) in this model does not obey cq. [S]on square lattice 
rather it's something like 



#(t)~logt. 



(6) 



This very weak chaos behavior, if one calls that so, is 
also confirmed on triangular (fig [3]) and rhombic lattices 
which shows that Manna's attractor is not alike BTW 
one, and they may not be classified in the same univer- 
sality class. 

Dhar-Ramaswamy model: figure shows the evolution 
of H(t) for this directed model. The few first points are 
attributed to lattice effects (as it's seen in other mod- 
els on different lattices too) and there is two different 
regimes of power behavior after which H(t) saturates be- 
cause of finite size effects; the exponents are 0.64 ± 0.04 
and 0.15 ±0.03 correspondingly. The intermediate regime 
is unknown to us. 

We have also studied a mixed Manna-Zhang model. In 



o 
u 

g 

'.7. 



60 
C 

6 

E 

X 




time 



FIG. 5: Evolution of H{t) for dissipative BTW model no 
square lattice with L = 256 for 5000 samples. At large times 
H(t) saturates because of finite size effects. The blue graph 
represents BTW without bulk dissipation. The green and red 
graphs correspond to 164(= 0.0025L 2 ) and 3277(= 0.05L 2 ) 
dissipative sites with fc = 1. The up-left graph shows how /3 
decreases rapidly by increase of dissipation amount (for fixed 
k(— 1)). The asymptotic value of /3 is 0.43. 
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FIG. 6: Evolution H(t) for dissipative Zhang model on square 
lattice for L — 256 and 10000 samples. The critical height is 
one and each site dissipates 0.5 unit of sand. In the main 
graph j3i = 0.64 and /?2 = 0.47. The up-left and down-right 
graphs show the rapid decrease of both /3i and ,82 by increase 
of dissipation amount. The asymptotic values are 0.62 and 
0.42 correspondingly 



this model when a site exceeds the critical height, cither 
the left and right or up and down neighbors are each given 
half of its energy by equal probability. Deviation from 
" weak chaos" behavior is also seen in this model. It seems 
that the stochastic property is responsible for transition 
from "weak chaos" to "very weak chaos" behavior. 

Since "weak chaos" is not observed in all sandpile mod- 
els, particularly in Manna model, self-organized criti- 
cality is not necessarily accompanied by "weak chaos". 
Therefore BTC's conjecture seems not to be correct in 
general. Since the "weak chaos" exponent of BTW and 
Zhang deterministic models do not depend on the lattice 
geometry, it may be considered as a universal property 
of these models and it may be viewed as a characteristic 
exponent (table HJ). If so, since BTW's weak chaos expo- 
nent is 0.74 ± 0.03 and Zhang's weak chaos exponent is 
0.96 ± 0.05 these two models (as an Abelian model and a 
non Abelian model) can't be classified in the same univer- 
sality class. In this manner our simulation based results 
do not agree with Zhang's conjecture about the unifica- 
tion of BTW model and his model in the thermodynamic 
limitii, Ben-Hur and Biham's classification^, and Daz- 
Guilera and Corral's classification^^. It should be noted 
that these classifications are all based on static consider- 
ations, where we have considered a dynamic property of 
the attractor which might cause this disagreement. 

On the other hand, although the root of this behavior 
of Manna model is not known, this classification based 
on the "weak chaos" exponent may put Manna model 
in a different class from BTW and Zhang models which 
is consistent from this aspect with Ben-Hur and Bi- 
ham's classification^ but inconsistent with Manna'p^ 1 
and Chcssa et al£- results. 



IV. DISSIPATIVE MODELS 

All wc have reported above is related to sandpile mod- 
els at the critical point. What happens to this behav- 
ior getting off the critical point? Do they still show 
weak chaos behavior? (It should be noted that since we 
have not studied fixed energy sandpile models, off critical 
states here means those states which their mean energy 
p satisfies p < p c , where p c is the critical mean energy) 

Adding bulk dissipation to sandpile models can be 
done using different methods. For discrete models the 
toppling rule is changed in some randomly chosen sites 
called dissipative sites. If the (i, j)th site is a dissipative 
one, it topples when z%j > z c + k where k is a posi- 
tive integer. In this method dissipation is controlled by 
two parameters, the number of dissipative sites and k. In 
continuous models there is no need of selected dissipative 
sites and we can impose the so called new toppling rule 
to all sites where k this time is a real number, therefore 
dissipation is controlled only by k. In discrete models wc 
can also use their continuous version and add arbitrary 
dissipation to all sites not to impose another stochastic 
parameter in the system. We have used the first method 
for the discrete BTW model and the second method for 
(continuous) Zhang model. 

Although bulk dissipation does not differ in principle 
from boundary dissipation, its importance is because it 
imposes a characteristic length in the system which de- 
stroys criticality (for example seeiS). Since bulk dissi- 
pation delays some toppling events we expect the rate of 
divergence of nearby configurations to decrease when dis- 
sipation increases. So we do not expect off critical BTW 
and Zhang models to show "weak chaos" behavior. 

Fig O shows that in BTW model by increase of dissi- 



pation, despite the interval of lattice effects which takes 
more long, the exponent of weak chaos region (3 decreases 
rapidly. This is exactly what we expect because the char- 
acteristic length rapidly becomes comparable with the 
system size, therefore "weak chaos" behavior decays. The 
up- left graph of figure [5] shows how rapid decrease of the 
exponent versus dissipation magnitude in BTW model. 

Fig [5] shows that Zhang model does not behave like 
BTW. The main graph shows the evolution of H(t) for 
0.5 sand unit dissipation at each site. Ascribing first few 
steps to lattice effects, there are to regimes of power- 
law behavior in this model which both exponents j3\ and 
P2 decrease rapidly as the magnitude of dissipation (the 
characteristic length of the system) increases. These are 
shown in up-left and down-right graphs of fig [6j Why 
two different regimes appear is steel unknown to us. 



stochastic property. Dhar-Ramaswamy directed model 
shows a more complicated behavior where it shows two 
independent regimes of power-law with an unknown 
regime in between. So our results contain both exam- 
ples and counter-examples of BTC's previous conjecture. 
Although weak chaos exponent does not seem to be a 
general characteristic of sandpile models but we have 
found a good evidence to concern it as a test for differ- 
ent universal classifications of these models. By means of 
this test (Abelian) BTW, (non-Abelian) Zhang, stochas- 
tic Manna, and directed Dhar-Ramaswamy models all 
seem to belong to different universality classes (where 
non of the offered classifications are in a complete accor- 
dance with). We have also shown that as we tend to off 
critical states, weak chaos behavior seems to disappear 
were in Zhang model an unknown split of H (t) into two 
power- law regimes is seen. 



V. CONCLUSIONS 

Our simulation results show that both BTW and 
Zhang models obey BTC's conjecture except some lat- 
tice effects where Manna model does not because of its 
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